F0. Pre-processing
Olga Viedma
2023-11-12
LIBRARIES
#1.CHM PITFREE 0.5 m
LIDAR_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
LIDAR_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (LIDAR_dir, pattern = "*.las", full.names = TRUE)
las_list1 <-list.files (LIDAR_dir, pattern = "*.las", full.names = FALSE)
short_name1<-stri_sub(las_list1, 1,-5)
lidar_maps<-lapply(las_list, function(X) lidR::readLAS(X, filter = "-drop_z_below 0"))
#plot(lidar_maps[[1]])
chm_pitfree_list<-list()
outputnames2<-NULL
for (i in seq_along(lidar_maps)) {
lidar_maps1<-lidar_maps[[i]]
chm_pitfree<- grid_canopy(lidar_maps1, res=0.5,pitfree( c(0,2,5,10,15,20,25,30,35,40), c(0,1.5), subcircle=0.15))
chm_pitfree[chm_pitfree > 40] <- NA
chm_pitfree[chm_pitfree < 0] <- 0
chm_pitfree1 <- projectRaster(chm_pitfree, crs=26916)
chm_pitfree_list[[i]]<-chm_pitfree1
outputnames2<- paste0(LIDAR_out ,"/", short_name1[i] ,"_CHM05m.tif",collapse = NULL,sep="")
writeRaster(chm_pitfree_list[[i]],outputnames2, overwrite=T)
}
col <- height.colors(25)
plot(chm_pitfree_list[[1]],col=col)
#2.DETECT TREE TOPS
ttop_dir<- file.path(system.file("extdata", package = "LadderFuelsR"))
las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*.las", full.names = TRUE)
las_list1 <-list.files (las_dir, pattern = "*.las", full.names = FALSE)
las_files<-lapply(las_list, function(X) lidR::readLAS(X, filter = "-drop_z_below 0"))
names1<-list()
filename_post2<-list()
chm_ttop_list2<-list()
for (j in seq_along(las_list)){
short_name<-stri_sub(las_list1[j], 1, -5)
names1<-rbind(names1,short_name)
filename_post2 <- paste0(names1,"_multi2_tree_tops" )
las_file <- lidR::readLAS(las_list[j],"-drop_z_below 0")
# parameters
ws= 2.5
hmin = 2
res=0.5
ttops_multichm = find_trees(las_file, multichm(res = res, dist_2d = 2,ws= ws, layer_thickness = 0.3,dist_3d = 1, hmin = hmin))
#ttops_multichm = find_trees(las_file, multichm(res = 2, dist_2d = 2,ws= 5, layer_thickness = 0.3,dist_3d = 1, hmin = 4))
proj4string(ttops_multichm) <- CRS('+init=EPSG:26916')
chm_ttop_list2[[j]]<-ttops_multichm
writeOGR(chm_ttop_list2[[j]], dsn=ttop_dir ,layer = filename_post2[[j]],driver = 'ESRI Shapefile',
overwrite_layer = T)
}
las1<-las_files[[1]]
chm1<-chm_pitfree_list[[1]]
ttops1<-chm_ttop_list2[[1]]
#plot(chm1, col = height.colors(50))
#plot(ttops1, add = TRUE, pch= 16, col = "black", main = "Tree tops over the CHM")
# Create an rgl point cloud
x<-add_treetops3d(plot(las_files[[1]], bg = "white", size = 4), ttops1)
# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 60, zoom = 0.75)
# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x", width = 800, height = 600)3D plot
#3.CROWNS (Silva) with locate_trees (ws: MULTICHM FIXED)
crowns_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
ttop_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
ttop_list <-list.files (ttop_dir, pattern = "*_multi2_tree_tops.shp", full.names = TRUE)
ttops_multichm<-lapply(ttop_list, function (x) st_read(x))## Reading layer `Eglin_zone1__crown2m_silva_multi2_tree_tops' from data source
## `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1__crown2m_silva_multi2_tree_tops.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 886 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 525912.8 ymin: 3378523 xmax: 526028.2 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
## Reading layer `Eglin_zone1_000000_multi2_tree_tops' from data source
## `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1_000000_multi2_tree_tops.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 886 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 525912.2 ymin: 3378523 xmax: 526028.2 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16Nchm_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
chm_list <-list.files (chm_dir, pattern = "_CHM05m.tif", full.names = TRUE)
las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*.las", full.names = TRUE)
names1<-list()
filename_post1<-list()
crowns_silva_list<-list()
for (j in seq_along(las_list)){
short_name<-gsub(".*/","",las_list[j])
short_name1<-stri_sub(short_name, 1, -11)
names1<-rbind(names1,short_name1)
filename_post1 <- paste0(names1,"crown2m_silva.las" )
chm05m <- raster(chm_list[j])
ttops<-ttops_multichm[[j]]
las_file <- lidR::readLAS(las_list[j])
algo_silva1 <-silva2016(chm05m, ttops, max_cr_factor = 0.6, exclusion = 0.3, ID = "treeID")
crowns_silva_las1 <-segment_trees(las_file, algo_silva1, attribute = "treeID", uniqueness = "incremental")
crowns_silva_las2<-filter_poi(crowns_silva_las1, !is.na(treeID))
crowns_silva_list[[j]]<-crowns_silva_las2
lidR::writeLAS(crowns_silva_list[[j]], file = paste0(crowns_dir, "/", filename_post1[[j]], sep=""))
}
my_palette <- colorRampPalette(col)
x1<-plot(crowns_silva_list[[1]], color = "treeID", pal = my_palette, bg = "white")
# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 10, zoom = 0.75)
# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x1", width = 800, height = 600)3D plot
FUNCTION CROWN METRICS
custom_crown_metrics <- function(z, i) { # user-defined function
metrics <- list(
dz = 1,
th = 1,
z_max = max(z), # max height
z_min = min(z), # min height
z_mean = mean(z), # mean height
z_sd = sd(z), # vertical variability of points
z_q1=quantile(z, probs = 0.01),
z_q5=quantile(z, probs = 0.05),
z_q25=quantile(z, probs = 0.25),
z_q50=quantile(z, probs = 0.50),
z_q75=quantile(z, probs = 0.75),
z_q95=quantile(z, probs = 0.95),
crr=(mean(z)-min(z))/(max(z)-min(z))
)
return(metrics) # output
}
ccm = ~custom_crown_metrics(z = Z, i = Intensity)#4.TREE METRICS CROWNS leafR (Silva)
stats_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
las_dir <- file.path(system.file("extdata", package = "leafR"))
las_list <-list.files (las_dir, pattern = "*_crown2m_silva.las", full.names = TRUE)
filename_post<-NULL
names2<-NULL
crown_list<-list()
for (j in seq_along(las_list)){
short_name<-stri_sub(las_list[j], 1, -5)
short_name1<-gsub(".*/","",short_name)
names2<-rbind(names2,short_name1)
filename_post <- paste0(names2,"_metrics_leafR.txt" )
normlas_file<-lidR::readLAS(las_list[[j]])
# Setting the xyz coordinates and subsetting the data
xyzi<-subset(normlas_file@data[,c(1:3,5,17)])
xyzi1<-filter(xyzi, Z >= 0.3)
xyzi1$treeID<-as.factor(xyzi1$treeID)
TreesMetrics<-CrownMetrics(xyzi1)
TreesMetrics1<-do.call(rbind, TreesMetrics)
TreesMetrics3 <- t(TreesMetrics1)
crown_list[[j]]<-TreesMetrics3
write.table(crown_list[[j]], paste0(stats_out, "/", filename_post[[j]] , sep=""), sep="\t", row.names = FALSE)
}## .................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................metrics1<-crown_list[[1]]
head(metrics1)## Tree TotalReturns ETOP NTOP EMIN NMIN EMAX NMAX
## [1,] 1 51 526023.5 3378613 526021.8 526024.9 3378611 3378614
## [2,] 2 136 525997.4 3378587 525995.8 526000 3378585 3378590
## [3,] 3 84 526020.6 3378529 526019.6 526022.4 3378527 3378532
## [4,] 4 100 525976.2 3378550 525974.9 525977.9 3378548 3378552
## [5,] 5 101 525965.9 3378556 525963 525967.6 3378554 3378557
## [6,] 6 325 525998.2 3378545 525997 526000.1 3378543 3378547
## EWIDTH NWIDTH HMAX HMIN HMEAN HMEDIAN HMODE HVAR HSD HCV HKUR HSKE
## [1,] 3.19 2.77 15.77 1.88 11.53 11.31 10.54 8.12 2.85 24.7 4.04 -0.71
## [2,] 4.12 4.66 14.99 3.13 11.28 11.47 14.51 7.85 2.8 24.83 2.35 -0.52
## [3,] 2.85 4.74 14.67 2.86 10.46 11.18 9.53 8.76 2.96 28.31 2.88 -0.88
## [4,] 2.98 4.32 14.66 0.57 10.1 8.84 5.4 13.69 3.7 36.62 2.41 -0.52
## [5,] 4.51 3.42 14.53 0.33 9.49 10.01 0.33 12.58 3.55 37.39 4.01 -1.04
## [6,] 3.09 3.98 14.53 0.33 10.36 11.62 13.62 12.67 3.56 34.34 3.88 -1.12
## H05TH H10TH H15TH H20TH H25TH H30TH H35TH H40TH H45TH H50TH H55TH H60TH
## [1,] 7.08 8.35 9.3 9.7 9.8 10.17 10.54 10.88 11.25 11.31 11.41 12.52
## [2,] 7.02 7.3 7.71 8.42 9.04 9.92 10.11 10.31 11.22 11.47 12.68 12.84
## [3,] 4.4 5.16 6.66 8.85 9.43 9.53 9.82 10.16 10.64 11.18 11.62 12.12
## [4,] 3.52 5.4 6.37 7.63 7.82 7.91 7.96 8.29 8.48 8.84 12.38 12.87
## [5,] 0.42 6.4 7.11 7.25 7.65 7.87 8.53 8.87 9.56 10.01 10.48 10.7
## [6,] 0.55 6.83 6.92 7.15 7.28 8.41 9.54 10.19 10.63 11.62 12.2 12.33
## H65TH H70TH H75TH H80TH H90TH H95TH H99TH IMAX IMIN IMEAN IMEDIAN IMODE
## [1,] 12.68 12.94 13.73 14.91 15.12 15.36 15.77 145 3 51.65 43 3
## [2,] 13.04 13.35 13.54 14.39 14.51 14.75 14.88 139 0 68.58 74.5 118
## [3,] 12.27 12.41 12.73 12.98 13.37 14.25 14.54 123 0 42.9 39.5 5
## [4,] 12.97 13.08 13.56 14 14.27 14.51 14.64 112 0 38.47 33 15
## [5,] 10.99 11.38 12.33 12.63 13.36 14.25 14.52 137 0 55.96 62 92
## [6,] 12.77 13.37 13.47 13.64 13.78 13.9 14.35 151 0 55.42 54 7
## IVAR ISD ICV IKUR ISKE I05TH I10TH I15TH I20TH I25TH I30TH I35TH
## [1,] 1464.79 38.27 74.1 1.95 0.35 3.5 7 8 12 17.5 22 28
## [2,] 1312.94 36.23 52.83 2.09 -0.28 6 12 21.5 33 43.75 49.5 55
## [3,] 1169.77 34.2 79.72 2.05 0.46 2 3.3 5 6.2 9.5 15.8 21.1
## [4,] 782.6 27.97 72.72 2.32 0.54 3.95 5.9 9.7 12.8 14.75 17 19
## [5,] 1387.94 37.26 66.57 1.8 0.06 2 6 11 14 21 29 35
## [6,] 1420.3 37.69 68.01 1.93 0.25 4 7 11 15 21 27 32
## I40TH I45TH I50TH I55TH I60TH I65TH I70TH I75TH I80TH I90TH I95TH I99TH
## [1,] 30 33 43 54 69 71.5 75 86 93.5 100 109.5 129
## [2,] 62 68.75 74.5 78 84 89 92 95.25 105.75 116 121.25 132.2
## [3,] 27.6 34.4 39.5 43 48.6 51.95 55 69.75 88 92.7 98.55 117.19
## [4,] 21 29.1 33 39 44.4 49 54.3 62.25 69.15 79.1 84.2 109.03
## [5,] 39 51 62 68 73 76 81 84 97 101 109 131
## [6,] 38 44.8 54 60 68 74 80 85 101 107.6 117.8 132.52#5.TREE METRICS CROWNS standard and customized (Silva)
stats_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*_crown2m_silva.las", full.names = TRUE)
names2<-list()
filename_post4<-list()
crowns_metrics_list<-list()
for (j in seq_along(las_list)){
short_name<-stri_sub(las_list[j], 1, -5)
short_name1<-gsub(".*/","",short_name)
names2<-rbind(names2,short_name1)
filename_post4 <- paste0(names2,"_metrics_convex.shp" )
las_file1 <- lidR::readLAS(las_list[j])
las_file2<-filter_poi(las_file1, !is.na(treeID))
las_file3<-filter_poi(las_file2, Z >= 1)
metrics1 = crown_metrics(las_file3,func = .stdtreemetrics, geom = "convex")
crown_diam<-data.frame(sqrt(metrics1$convhull_area/ pi) * 2)
names(crown_diam)<-"crown_diam"
metrics2 = crown_metrics(las_file3,func = ccm, geom = "convex") #concave
metrics_all <- dplyr::bind_cols(list(metrics1,crown_diam,metrics2))
metrics_all1 <- metrics_all[,c(1:4,6,10:21)]
names(metrics_all1)<-c("treeID", "Z", "npoints", "convhull_area", "crown_diam", "z_max", "z_min", "z_mean","z_sd", "z_q1","z_q5", "z_q25","z_q50","z_q75", "z_q95", "crr", "geometry" )
metrics_all2<- st_as_sf(metrics_all1)
crowns_metrics_list[[j]]<-metrics_all2
st_write(crowns_metrics_list[[j]], paste0(dsn=stats_out, "/", filename_post4[[j]],sep=""), driver = 'ESRI Shapefile',append=F)
}## Writing layer `Eglin_zone1__crown2m_crown2m_silva_metrics_convex' to data source
## `D:/OLGA/R_library/LadderFuelsR/extdata/Eglin_zone1__crown2m_crown2m_silva_metrics_convex.shp' using driver `ESRI Shapefile'
## Writing 878 features with 16 fields and geometry type Polygon.
## Deleting layer `Eglin_zone1__crown2m_silva_metrics_convex' using driver `ESRI Shapefile'
## Writing layer `Eglin_zone1__crown2m_silva_metrics_convex' to data source
## `D:/OLGA/R_library/LadderFuelsR/extdata/Eglin_zone1__crown2m_silva_metrics_convex.shp' using driver `ESRI Shapefile'
## Writing 877 features with 16 fields and geometry type Polygon.ttops1<-st_as_sf(chm_ttop_list2[[1]])
crowns1<-st_as_sf(crowns_metrics_list[[1]])
ttops_within_crowns <- st_intersection(ttops1, crowns1)
# Set the size of the plotting device
par(mfrow = c(1, 1), mar = c(1, 1, 1, 1), pin = c(5, 4))
plot(st_geometry(crowns1), pch = 16, col = "green")
plot(ttops_within_crowns, add = TRUE, pch= 16, col = "darkblue", main = "Tree tops over the crowns")
#6.NO OVERLAPPING CROWN POLYGONS
dir_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
tree_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
tree_list <-list.files (tree_dir, pattern = "*_crown2m_silva_metrics_convex.shp", full.names = TRUE)
tree_files <-lapply(tree_list, function (X) st_read(X))## Reading layer `Eglin_zone1__crown2m_crown2m_silva_metrics_convex' from data source `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1__crown2m_crown2m_silva_metrics_convex.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 878 features and 16 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 525912.5 ymin: 3378524 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N
## Reading layer `Eglin_zone1__crown2m_silva_metrics_convex' from data source
## `D:\OLGA\R_library\LadderFuelsR\extdata\Eglin_zone1__crown2m_silva_metrics_convex.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 877 features and 16 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 525912.5 ymin: 3378523 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16Npolys1<-tree_files[[1]]
for(i in 1:nrow(polys1)) {
other_polys <- st_union(polys1[-i, ]$geometry)
polys1[i, ]$geometry <- st_difference(polys1[i, ]$geometry, other_polys)
}
st_write(polys1,paste0(tree_dir, "/", "crown2m_silva_convex_no_overlap.shp"), append=FALSE)## Deleting layer `crown2m_silva_convex_no_overlap' using driver `ESRI Shapefile'
## Writing layer `crown2m_silva_convex_no_overlap' to data source
## `D:/OLGA/R_library/LadderFuelsR/extdata/crown2m_silva_convex_no_overlap.shp' using driver `ESRI Shapefile'
## Writing 878 features with 16 fields and geometry type Polygon.# Set the size of the plotting device
par(mfrow = c(1, 1), mar = c(1, 1, 1, 1), pin = c(5, 4))
plot(st_geometry(polys1), pch = 16, col = "orange")
plot(ttops_within_crowns, add = TRUE, pch= 16, col = "blue", main = "Tree tops over the crowns")
#7.CROP LAS FILES WITH CROWN POLYGONS
dir_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
las_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
las_list <-list.files (las_dir, pattern = "*_crown2m_silva.las", full.names = TRUE)
las_files <-lapply(las_list, function (X) lidR::readLAS(X))
tree_dir <- file.path(system.file("extdata", package = "LadderFuelsR"))
tree_list <-list.files (tree_dir, pattern = "*_no_overlap.shp", full.names = TRUE)
tree_files <-lapply(tree_list, function (X) st_read(X))## Reading layer `crown2m_silva_convex_no_overlap' from data source
## `D:\OLGA\R_library\LadderFuelsR\extdata\crown2m_silva_convex_no_overlap.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 878 features and 16 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 525912.5 ymin: 3378524 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16N###################################
names1 <- list()
filename_post3 <- list()
for (j in seq_along(las_list)) {
short_name <- stri_sub(las_list[j])
short_name1 <- gsub(".*/", "", short_name)
short_name2 <- stri_sub(short_name1, 1, 11)
print(short_name2)
if (is.na(short_name2)) {
# Handle NA values or skip them
next
}
names1[[j]] <- short_name2
filename_post3[[j]] <- paste0(short_name2, "_CROWN.las")
las1 <- las_files[[j]]
trees1 <- tree_files[[j]]
trees_ID <- trees1 %>% dplyr::select(treeID)
n <- nrow(trees_ID)
crown_cort <- vector("list", length=n)
for (i in 1:n) {
kk <- trees_ID[i,]
crown_cort[[i]] = clip_roi(las1, kk)
lidR::writeLAS(crown_cort[[i]], paste0(dir_out,"/", i, "_", filename_post3[j]))
}
}## [1] "Eglin_zone1"my_palette <- colorRampPalette(col)
x2<-plot(crown_cort[[1]], color = "Z", pal = my_palette, bg = "black", size = 2.5)
# Customize the plot orientation
rgl.viewpoint(theta = 0, phi = 0, fov = 60, zoom = 0.75)
# Convert the rgl scene to an HTML widget
rglwidget(elementId = "x2", width = 400, height = 600)3D plot
#8.LAI-LAD METRICS BY TREE
stats_out<-file.path(system.file("extdata", package = "LadderFuelsR"))
las_dir1 <- file.path(system.file("extdata", package = "leafR"))
las_list1 <- list.files(las_dir1, pattern = "*_CROWN.las", full.names = TRUE, ignore.case = TRUE)
# create a vector to hold the file names of .las files with more than 3 points
files_with_more_than_3_points <- c()
# loop through each file
for (file in las_list1) {
las_data <- lidR::readLAS(file)
las_data1<-filter_poi(las_data, Z >= 1)
# skip to next file if there was a problem reading
if (is.null(las_data1)) next
# check if it contains more than three points
if (las_data1@header$`Number of point records` > 3) {
files_with_more_than_3_points <- c(files_with_more_than_3_points, file)
}
}
# Creates a data frame of the 3D voxels information (xyz) with Leaf Area Density values
short_name1<-NULL
profile_list<-NULL
lidar_lai_list<-NULL
understory_lai_list<-NULL
LAHV_metric_list<-NULL
for (j in seq_along(files_with_more_than_3_points)){
short_name<-stri_sub(files_with_more_than_3_points[j], 1, -5)
short_name1<-gsub(".*/","",short_name)
normlas_file<-files_with_more_than_3_points[[j]]
VOXELS_LAD = lad.voxels(normlas_file, grain.size = 2)
lad_profile = lad.profile(VOXELS_LAD, relative = F)
lai_tot = lai(lad_profile)
understory_lai <- lai(lad_profile, min = 0.3, max = 2.5)
LAHV_metric<- LAHV(lad_profile, LAI.weighting = FALSE, height.weighting = FALSE) # Leaf Area height volume (wood volume, AGB)
lad_profile1 = data.frame(lad_profile, treeID = short_name1)
lai_tot1 = data.frame(lai_tot, treeID = short_name1)
understory_lai1 = data.frame(understory_lai, treeID = short_name1)
LAHV_metric1 = data.frame(LAHV_metric, treeID = short_name1)
profile_list<-rbind(profile_list, lad_profile1)
lidar_lai_list<-rbind(lidar_lai_list,lai_tot1)
understory_lai_list <-rbind(understory_lai_list,understory_lai1)
LAHV_metric_list<-rbind(LAHV_metric_list,LAHV_metric1)
}
write.table(profile_list,file =paste(stats_out, "/","alltrees_LAD_profile_voxels2m",".txt", sep=""), sep="\t", row.names = FALSE)
write.table(lidar_lai_list,file =paste(stats_out, "/","alltrees_LAI_profile_voxels2m_",".txt", sep=""), sep="\t", row.names = FALSE)
write.table(understory_lai_list,file =paste(stats_out, "/","alltrees_understory_profile_voxels2m",".txt", sep=""), sep="\t", row.names = FALSE)
write.table(LAHV_metric_list,file =paste(stats_out, "/","alltrees_LAHV_profile_voxels2m",".txt", sep=""), sep="\t", row.names = FALSE)
head(profile_list)## height lad treeID
## 1 1.5 0.01724822 1_Eglin_zone1_CROWN
## 2 2.5 0.00000000 1_Eglin_zone1_CROWN
## 3 3.5 0.00000000 1_Eglin_zone1_CROWN
## 4 4.5 0.00000000 1_Eglin_zone1_CROWN
## 5 5.5 0.00000000 1_Eglin_zone1_CROWN
## 6 6.5 0.10136628 1_Eglin_zone1_CROWN#9.DEPURATING TREE LAD PROFILES (> 3 HEIGHT VALUES)
LAD_path<- file.path(system.file("extdata", package = "LadderFuelsR"),"alltrees_LAD_profile_voxels2m.txt")
stats_tot<- read.table(LAD_path,sep="\t", header=T)
cols <- c('treeID')
stats_tot[cols] <- lapply(stats_tot[cols], function (x) as.factor(x))
stats_tot$lad<-round(stats_tot$lad,digits = 4)
cases <- data.frame(table(stats_tot$treeID))
cases1 <-cases[cases$Freq > 3, ]
names(cases1)<-c("treeID", "Freq")
stats_tot1 <- stats_tot[stats_tot$treeID %in% cases1$treeID, ]
stats_tot2 <- data.frame(stats_tot1 %>% replace(is.na(.), 0.01))
output_file <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
write.table(stats_tot2, file = output_file, sep = "\t", row.names = FALSE)#10.GAPS AND FUEL LAYERS BASE HEIGHT (FBH)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)
trees_name1 <- as.character(LAD_profiles$treeID)
trees_name2 <- factor(unique(trees_name1))
metrics_precentile_list<-list()
for (i in levels(trees_name2)) {
tree2 <- LAD_profiles |> dplyr::filter(treeID == i)
metrics_precentil <- get_gaps_fbhs(tree2)
metrics_precentile_list[[i]] <- metrics_precentil
}
metrics_all_percentil <- dplyr::bind_rows(metrics_precentile_list)
metrics_all_percentil$treeID <- factor(metrics_all_percentil$treeID)
# Remove the row with all NA values from the original data frame
# First remove "treeID" and "treeID1" columns
tree_metrics_no_treeID <- metrics_all_percentil[, -which(names(metrics_all_percentil) == c("treeID","treeID1"))]
# Check if any row has all NA values
rows_with_all_NA_or_zero <- apply(tree_metrics_no_treeID, 1, function(row) all(is.na(row) | row == 0))
# Get the row index with all NA values
row_index <- which(rows_with_all_NA_or_zero)
# Remove the row with all NA values from the original data frame
if (length(row_index) > 0) {
tree_metrics_filtered <- metrics_all_percentil[-row_index, ]
} else {
tree_metrics_filtered <- metrics_all_percentil
}
write.table(tree_metrics_filtered, file= file.path(system.file("extdata", package = "LadderFuelsR"), "1_gaps_fbhs_metrics.txt"),
sep = "\t",row.names = FALSE)
head(tree_metrics_filtered)## treeID cbh1 gap1 gap2 cbh2 cbh3 cbh4 gap_lad1
## height...1 1_Eglin_zone1_CROWN 1.5 2.5 5.5 6.5 9.5 14.5 0
## height...2 10_Eglin_zone1_CROWN 5.5 1.5 4.5 10.5 13.5 <NA> 0
## height...3 100_Eglin_zone1_CROWN 1.5 2.5 4.5 5.5 8.5 11.5 0
## height...4 101_Eglin_zone1_CROWN 4.5 1.5 3.5 6.5 11.5 <NA> 0
## height...5 102_Eglin_zone1_CROWN 5.5 1.5 4.5 8.5 11.5 <NA> 0
## height...6 103_Eglin_zone1_CROWN 1.5 2.5 4.5 5.5 7.5 10.5 0
## gap_lad2 cbh_perc1 cbh_perc2 cbh_perc3 cbh_perc4 cbh_lad1 cbh_lad2
## height...1 0 50 75 90 50 0.0172 0.1014
## height...2 0 75 75 100 NA 0.2310 0.0608
## height...3 0 50 50 95 75 0.0239 0.0143
## height...4 0 75 100 50 NA 0.1032 0.2544
## height...5 0 90 100 50 NA 0.1754 0.2338
## height...6 0 50 75 95 100 0.0230 0.0445
## cbh_lad3 cbh_lad4 max_height treeID1 gap3 gap_lad3 cbh5 cbh_perc5
## height...1 0.1933 0.0333 15.5 1 <NA> NA <NA> NA
## height...2 0.6150 NA 14.5 10 9.5 0 <NA> NA
## height...3 0.1705 0.0924 12.5 100 <NA> NA <NA> NA
## height...4 0.0668 NA 12.5 101 <NA> NA <NA> NA
## height...5 0.0274 NA 12.5 102 <NA> NA <NA> NA
## height...6 0.1744 0.3261 12.5 103 <NA> NA 11.5 50
## cbh_lad5 cbh6 cbh_perc6 cbh_lad6 gap4 gap5 gap_lad4 gap_lad5
## height...1 NA <NA> NA NA <NA> <NA> NA NA
## height...2 NA <NA> NA NA <NA> <NA> NA NA
## height...3 NA <NA> NA NA <NA> <NA> NA NA
## height...4 NA <NA> NA NA <NA> <NA> NA NA
## height...5 NA <NA> NA NA <NA> <NA> NA NA
## height...6 0.0343 <NA> NA NA <NA> <NA> NA NA#11.DISTANCE BETWEEN FUEL LAYERS
# Tree metrics derived from get_gaps_fbhs() function
gaps_fbhs_metrics_path <- system.file("extdata", "1_gaps_fbhs_metrics.txt", package = "LadderFuelsR")
gaps_fbhs_metrics <- read.table(gaps_fbhs_metrics_path, sep="\t", header=TRUE)
gaps_fbhs_metrics$treeID <- factor(gaps_fbhs_metrics$treeID)
trees_name1 <- as.character(gaps_fbhs_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))
metrics_distance_list <- list()
for (i in levels(trees_name2)) {
# Filter data for each tree
tree2 <- gaps_fbhs_metrics |> dplyr::filter(treeID == i)
# Get distance metrics for each tree
metrics_distance <- get_distance(tree2)
metrics_distance_list[[i]] <- metrics_distance
}
# Combine the individual data frames
metrics_all_distance <- dplyr::bind_rows(metrics_distance_list)
distance_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "2_distance_metrics.txt")
write.table(metrics_all_distance, file=distance_path, sep="\t", row.names=FALSE)
metrics_all_distance <- metrics_all_distance[, order(names(metrics_all_distance))]
head(metrics_all_distance)## cbh1 cbh2 cbh3 cbh4 cbh5 cbh6 dist1 dist2 dist3 gap1 gap2 gap3 gap4 gap5
## 1 1.5 6.5 9.5 14.5 NA NA 4 NA NA 2.5 5.5 NA NA NA
## 2 5.5 10.5 13.5 NA NA NA 4 1 NA 1.5 4.5 9.5 NA NA
## 3 1.5 5.5 8.5 11.5 NA NA 3 NA NA 2.5 4.5 NA NA NA
## 4 4.5 6.5 11.5 NA NA NA 3 NA NA 1.5 3.5 NA NA NA
## 5 5.5 8.5 11.5 NA NA NA 4 NA NA 1.5 4.5 NA NA NA
## 6 1.5 5.5 7.5 10.5 11.5 NA 3 NA NA 2.5 4.5 NA NA NA
## Hdist1 Hdist2 Hdist3 max_height treeID treeID1
## 1 5.5 NA NA 15.5 1_Eglin_zone1_CROWN 1
## 2 4.5 9.5 NA 14.5 10_Eglin_zone1_CROWN 10
## 3 4.5 NA NA 12.5 100_Eglin_zone1_CROWN 100
## 4 3.5 NA NA 12.5 101_Eglin_zone1_CROWN 101
## 5 4.5 NA NA 12.5 102_Eglin_zone1_CROWN 102
## 6 4.5 NA NA 12.5 103_Eglin_zone1_CROWN 103#12.FUEL LAYERS DEPTH
library(dplyr)
library(magrittr)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)
# Tree metrics derived from get_distance() function
distance_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "2_distance_metrics.txt")
distance_metrics <- read.table(distance_path, sep = "\t", header = TRUE)
distance_metrics$treeID <- factor(distance_metrics$treeID)
metrics_depth_list <- list()
for (i in levels(LAD_profiles$treeID)){
tree1 <- LAD_profiles |> dplyr::filter(treeID == i)
tree2 <- distance_metrics |> dplyr::filter(treeID == i)
# Get depths for each tree
metrics_depth <- get_depths(tree1, tree2)
metrics_depth_list[[i]] <- metrics_depth
}
# Combine the individual data frames
metrics_all_depth <- dplyr::bind_rows(metrics_depth_list)
depth_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "3_depth_metrics.txt")
write.table(metrics_all_depth, file = depth_path, sep = "\t", row.names = FALSE)
metrics_all_depth <- metrics_all_depth[, order(names(metrics_all_depth))]
head(metrics_all_depth)## cbh1 cbh2 cbh3 cbh4 cbh5 cbh6 depth0 depth1 depth2 depth3 dist1 dist2 dist3
## 1 1.5 6.5 9.5 14.5 NA NA 1 8 NA NA 4 NA NA
## 2 5.5 10.5 13.5 NA NA NA 0 4 3 NA 4 1 NA
## 3 1.5 5.5 8.5 11.5 NA NA 1 6 NA NA 3 NA NA
## 4 4.5 6.5 11.5 NA NA NA 0 7 NA NA 3 NA NA
## 5 5.5 8.5 11.5 NA NA NA 0 6 NA NA 4 NA NA
## 6 1.5 5.5 7.5 10.5 11.5 NA 1 6 NA NA 3 NA NA
## gap1 gap2 gap3 gap4 gap5 Hdepth0 Hdepth1 Hdepth2 Hdepth3 Hdist1 Hdist2 Hdist3
## 1 2.5 5.5 NA NA NA 1.5 14.5 NA NA 5.5 NA NA
## 2 1.5 4.5 9.5 NA NA 0.5 8.5 13.5 NA 4.5 9.5 NA
## 3 2.5 4.5 NA NA NA 1.5 11.5 NA NA 4.5 NA NA
## 4 1.5 3.5 NA NA NA 0.5 11.5 NA NA 3.5 NA NA
## 5 1.5 4.5 NA NA NA 0.5 11.5 NA NA 4.5 NA NA
## 6 2.5 4.5 NA NA NA 1.5 11.5 NA NA 4.5 NA NA
## max_height treeID treeID1
## 1 15.5 1_Eglin_zone1_CROWN 1
## 2 14.5 10_Eglin_zone1_CROWN 10
## 3 12.5 100_Eglin_zone1_CROWN 100
## 4 12.5 101_Eglin_zone1_CROWN 101
## 5 12.5 102_Eglin_zone1_CROWN 102
## 6 12.5 103_Eglin_zone1_CROWN 103#13.PLOTS GAPS AND FUEL LAYERS BASE HEIGHT (FBH)
library(LadderFuelsR)
library(ggplot2)
library(lattice)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)
# Tree metrics derived from get_depths() function
depth_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "3_depth_metrics.txt")
depth_metrics <- read.table(depth_path, sep = "\t", header = TRUE)
depth_metrics$treeID <- factor(depth_metrics$treeID)
# Generate plots for gaps and fbhs
plots_gaps_fbhs <- get_plots_gap_fbh(LAD_profiles, depth_metrics)
# Save plots for each tree
for (name in names(plots_gaps_fbhs)) {
#print(plots_gaps_fbhs[[name]]) # print the plot
ggsave(file.path(system.file("extdata", package = "LadderFuelsR"), paste0( name, "_gap_fbh", ".tiff")), plot = plots_gaps_fbhs[[name]])
}
par(mfrow = c(3, 1))
# Plot the first panel
plot(plots_gaps_fbhs[[1]],width = 5, height = 5)
# Plot the second panel
plot(plots_gaps_fbhs[[2]],width = 5, height = 5)
# Plot the third panel
plot(plots_gaps_fbhs[[3]],width = 5, height = 5)
#14.FUEL LAYERS BASE HEIGHT (FBH) AFTER REMOVING DISTANCES = 1
library(SSBtools)
library(dplyr)
library(magrittr)
# Tree metrics derived from get_depths() function
depth_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "3_depth_metrics.txt")
depth_metrics <- read.table(depth_path, sep = "\t", header = TRUE)
depth_metrics$treeID <- factor(depth_metrics$treeID)
trees_name1 <- as.character(depth_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))
fbh_corr_list <- list()
for (i in levels(trees_name2)){
# Filter data for each tree
tree3 <- depth_metrics |> dplyr::filter(treeID == i)
# Get real fbh for each tree
fbh_corr <- get_real_fbh(tree3)
# Store fbh values in a list
fbh_corr_list[[i]] <- fbh_corr
}
# Combine fbh values for all trees
fbh_corr_all <- dplyr::bind_rows(fbh_corr_list)
fbh_corr_all$treeID <- factor(fbh_corr_all$treeID)
# Reorder columns
# Get original column names
original_column_names <- colnames(fbh_corr_all)
# Specify prefixes
prefixes <- c("treeID", "Hdist", "Hcbh", "Hdepth", "dist", "depth", "max_height")
# Initialize vector to store new order
new_order <- c()
# Loop over prefixes
for (prefix in prefixes) {
# Find column names matching the current prefix
matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)
# Append to the new order
new_order <- c(new_order, matching_columns)
}
# Reorder values
fbh_corr_all <- fbh_corr_all[, new_order]
# Save the reordered data
fbh_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "4_fbh_metrics_corr.txt")
write.table(fbh_corr_all, file = fbh_corr_path, sep = "\t", row.names = FALSE)
head(fbh_corr_all)## treeID treeID1 Hdist1 Hdist2 Hdist3 Hcbh1 Hcbh2
## new_Hcbh...1 1_Eglin_zone1_CROWN 1 5.5 NA NA 1.5 6.5
## new_Hcbh...2 10_Eglin_zone1_CROWN 10 4.5 9.5 NA 5.5 5.5
## new_Hcbh...3 100_Eglin_zone1_CROWN 100 4.5 NA NA 1.5 5.5
## new_Hcbh...4 101_Eglin_zone1_CROWN 101 3.5 NA NA 4.5 4.5
## new_Hcbh...5 102_Eglin_zone1_CROWN 102 4.5 NA NA 5.5 5.5
## new_Hcbh...6 103_Eglin_zone1_CROWN 103 4.5 NA NA 1.5 5.5
## Hcbh3 Hcbh4 Hcbh5 Hcbh6 Hdepth1 Hdepth2 Hdepth3 Hdepth4 dist1
## new_Hcbh...1 6.5 6.5 NA NA 1.5 14.5 NA NA 4
## new_Hcbh...2 5.5 NA NA NA 8.5 13.5 NA NA 4
## new_Hcbh...3 5.5 5.5 NA NA 1.5 11.5 NA NA 3
## new_Hcbh...4 4.5 NA NA NA 11.5 NA NA NA 3
## new_Hcbh...5 5.5 NA NA NA 11.5 NA NA NA 4
## new_Hcbh...6 5.5 5.5 5.5 NA 1.5 11.5 NA NA 3
## dist2 dist3 depth1 depth2 depth3 depth4 max_height
## new_Hcbh...1 NA NA 1 8 NA NA 15.5
## new_Hcbh...2 1 NA 4 3 NA NA 14.5
## new_Hcbh...3 NA NA 1 6 NA NA 12.5
## new_Hcbh...4 NA NA 7 NA NA NA 12.5
## new_Hcbh...5 NA NA 6 NA NA NA 12.5
## new_Hcbh...6 NA NA 1 6 NA NA 12.5#15.FUEL LAYERS DEPTH AFTER REMOVING DISTANCES = 1
library(dplyr)
library(magrittr)
library(tidyr)
# Tree metrics derived from get_real_fbh() function
fbhcor_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "4_fbh_metrics_corr.txt")
fbh_metrics_corr <- read.table(fbhcor_path, sep = "\t", header = TRUE)
fbh_metrics_corr$treeID <- factor(fbh_metrics_corr$treeID)
trees_name1 <- as.character(fbh_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))
depth_metrics_corr <- lapply(levels(trees_name2), function(i) {
# Filter data for each tree
tree2 <- fbh_metrics_corr |> dplyr::filter(treeID == i)
# Get real depths for each tree
get_real_depths(tree2)
})
# Combine depth values for all trees
depth_corr_all <- dplyr::bind_rows(depth_metrics_corr)
depth_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "5_tree_depth_metrics_corr.txt")
write.table(depth_corr_all, file = depth_corr_path, sep = "\t", row.names = FALSE)
head(depth_corr_all)## treeID1 treeID Hcbh1 Hcbh2 dptf1 dptf2 Hdptf1 Hdptf2 dist1
## 1 1 1_Eglin_zone1_CROWN 1.5 6.5 1 8 1.5 14.5 4
## 2 10 10_Eglin_zone1_CROWN 5.5 NA 8 NA 13.5 NA 4
## 3 100 100_Eglin_zone1_CROWN 1.5 5.5 1 6 1.5 11.5 3
## 4 101 101_Eglin_zone1_CROWN 4.5 NA 7 NA 11.5 NA 3
## 5 102 102_Eglin_zone1_CROWN 5.5 NA 6 NA 11.5 NA 4
## 6 103 103_Eglin_zone1_CROWN 1.5 5.5 1 6 1.5 11.5 3
## Hdist1 max_height Hcbh3 dptf3 Hdptf3 dist2 dist3 Hdist2 Hdist3
## 1 5.5 15.5 NA NA NA NA NA NA NA
## 2 4.5 14.5 NA NA NA NA NA NA NA
## 3 4.5 12.5 NA NA NA NA NA NA NA
## 4 3.5 12.5 NA NA NA NA NA NA NA
## 5 4.5 12.5 NA NA NA NA NA NA NA
## 6 4.5 12.5 NA NA NA NA NA NA NA#16.FUEL LAYERS DISTANCES (> 1 M) AND CBH BASE ON MAXIMUM- AND LAST- DISTANCE
library(dplyr)
library(magrittr)
library(stringr)
# Tree metrics derived from get_real_depths() function
depth_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "5_tree_depth_metrics_corr.txt")
depth_metrics_corr <- read.table(depth_corr_path, sep = "\t", header = TRUE)
depth_metrics_corr$treeID <- factor(depth_metrics_corr$treeID)
trees_name1 <- as.character(depth_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))
distance_metrics_corr <- lapply(levels(trees_name2), function(i) {
# Filter data for each tree
tree2 <- depth_metrics_corr |> dplyr::filter(treeID == i)
# Get effective gap for each tree
get_effective_gap(tree2)
})
# Combine the individual data frames
distances_corr_all <- dplyr::bind_rows(distance_metrics_corr)
# =======================================================================#
# REORDER COLUMNS:
# =======================================================================#
# Get original column names
original_column_names <- colnames(distances_corr_all)
# Specify prefixes
prefixes <- c("treeID", "Hcbh", "dptf", "Hdptf", "effdist", "dist", "Hdist", "max_Hcbh", "max_dptf", "max_Hdptf", "last_Hcbh", "last_dptf", "last_Hdptf", "max_height")
# Initialize vector to store new order
new_order <- c()
# Loop over prefixes
for (prefix in prefixes) {
# Find column names matching the current prefix
matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)
# Extract numeric suffixes and order the columns based on these suffixes
numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
matching_columns <- matching_columns[order(numeric_suffixes)]
# Append to new order
new_order <- c(new_order, matching_columns)
}
# Reorder values
distances_corr_all1 <- distances_corr_all[, new_order]
distance_corr_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "6_tree_distances_metrics_corr.txt")
write.table(distances_corr_all1, file = distance_corr_path, sep = "\t", row.names = FALSE)
head(distances_corr_all1)## treeID1 treeID Hcbh1 Hcbh2 Hcbh3 dptf1 dptf2 dptf3 Hdptf1
## 1 1 1_Eglin_zone1_CROWN 1.5 6.5 NA 1 8 NA 1.5
## 2 10 10_Eglin_zone1_CROWN 5.5 NA NA 8 NA NA 13.5
## 3 100 100_Eglin_zone1_CROWN 1.5 5.5 NA 1 6 NA 1.5
## 4 101 101_Eglin_zone1_CROWN 4.5 NA NA 7 NA NA 11.5
## 5 102 102_Eglin_zone1_CROWN 5.5 NA NA 6 NA NA 11.5
## 6 103 103_Eglin_zone1_CROWN 1.5 5.5 NA 1 6 NA 1.5
## Hdptf2 Hdptf3 effdist1 effdist2 dist1 dist2 dist3 Hdist1 Hdist2 Hdist3
## 1 14.5 NA 4 NA 4 NA NA 5.5 NA NA
## 2 NA NA 4 NA 4 NA NA 4.5 NA NA
## 3 11.5 NA 3 NA 3 NA NA 4.5 NA NA
## 4 NA NA 3 NA 3 NA NA 3.5 NA NA
## 5 NA NA 4 NA 4 NA NA 4.5 NA NA
## 6 11.5 NA 3 NA 3 NA NA 4.5 NA NA
## max_Hcbh max_dptf max_Hdptf last_Hcbh last_dptf last_Hdptf max_height
## 1 6.5 8 14.5 6.5 8 14.5 15.5
## 2 5.5 8 13.5 5.5 8 13.5 14.5
## 3 5.5 6 11.5 5.5 6 11.5 12.5
## 4 4.5 7 11.5 4.5 7 11.5 12.5
## 5 5.5 6 11.5 5.5 6 11.5 12.5
## 6 5.5 6 11.5 5.5 6 11.5 12.5#17.FUELS LAD PERCENTAGE (> 25 %) AND CBH BASED ON MAXIMUM LAD PERCENTAGE
library(dplyr)
library(magrittr)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- system.file("extdata", "LAD_profiles.txt", package = "LadderFuelsR")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)
# Tree metrics derived from get_effective_gap() function
distcor_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "6_tree_distances_metrics_corr.txt")
distances_metrics_corr <- read.table(distcor_path, sep = "\t", header = TRUE)
distances_metrics_corr$treeID <- factor(distances_metrics_corr$treeID)
trees_name1 <- as.character(distances_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))
LAD_metrics1 <- list()
LAD_metrics2 <- list()
for (i in levels(trees_name2)) {
# Filter data for each tree
tree1 <- LAD_profiles |> dplyr::filter(treeID == i)
tree2 <- distances_metrics_corr |> dplyr::filter(treeID == i)
# Get LAD metrics for each tree
LAD_metrics <- get_layers_lad(tree1, tree2)
LAD_metrics1[[i]] <- LAD_metrics$df1
LAD_metrics2[[i]] <- LAD_metrics$df2
}
LAD_metrics_all1 <- dplyr::bind_rows(LAD_metrics1)
LAD_metrics_all2 <- dplyr::bind_rows(LAD_metrics2)
# List of data frames
LAD_metrics_list <- list(LAD_metrics_all1, LAD_metrics_all2)
# Initialize an empty list to store reordered data frames
reordered_list <- list()
# Specify prefixes (adjust accordingly)
prefixes <- c("treeID", "Hdist", "Hcbh", "effdist", "dptf", "Hdptf", "max", "last")
# Loop over each data frame
for (i in seq_along(LAD_metrics_list)) {
LAD_metrics_all <- LAD_metrics_list[[i]]
# Get original column names
original_column_names <- colnames(LAD_metrics_all)
# Initialize vector to store new order
new_order <- c()
# Loop over prefixes
for (prefix in prefixes) {
# Find column names matching the current prefix
matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)
# Extract numeric suffixes and order the columns based on these suffixes
numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
# Order the columns based on numeric suffixes
matching_columns <- matching_columns[order(numeric_suffixes)]
# Append to new order
new_order <- c(new_order, matching_columns)
}
# Reorder columns
LAD_metrics_all <- LAD_metrics_all[, new_order]
# Store the reordered data frame in the list
reordered_list[[i]] <- LAD_metrics_all
}
# Write the reordered data frame to a file
fuels_lad_all_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_all.txt")
write.table(reordered_list[[1]], file = fuels_lad_all_path, sep = "\t", row.names = FALSE)
fuels_lad_gt25_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_gt25perc.txt")
write.table(reordered_list[[2]], file = fuels_lad_gt25_path, sep = "\t", row.names = FALSE)
head(reordered_list[[2]])## treeID1 treeID Hdist1 Hdist2 Hcbh1 Hcbh2
## Percentage...1 1 1_Eglin_zone1_CROWN 5.5 NA 6.5 NA
## Percentage...2 10 10_Eglin_zone1_CROWN 4.5 NA 5.5 NA
## Percentage...3 100 100_Eglin_zone1_CROWN 4.5 NA 5.5 NA
## Percentage...4 101 101_Eglin_zone1_CROWN 3.5 NA 4.5 NA
## Percentage...5 102 102_Eglin_zone1_CROWN 4.5 NA 5.5 NA
## Percentage...6 103 103_Eglin_zone1_CROWN 4.5 NA 5.5 NA
## Hcbh1_Hdptf1 Hcbh2_Hdptf2 effdist1 dptf1 dptf2 Hdptf1 Hdptf2
## Percentage...1 94.49448 NA 5 8 NA 14.5 NA
## Percentage...2 99.37928 NA 4 8 NA 13.5 NA
## Percentage...3 94.88224 NA 4 6 NA 11.5 NA
## Percentage...4 99.66928 NA 3 7 NA 11.5 NA
## Percentage...5 98.91892 NA 4 6 NA 11.5 NA
## Percentage...6 95.17946 NA 4 6 NA 11.5 NA
## max_Hcbh max_dptf max_Hdptf max_height maxlad_Hcbh maxlad_Hdist
## Percentage...1 6.5 8 14.5 15.5 6.5 5.5
## Percentage...2 5.5 8 13.5 14.5 5.5 4.5
## Percentage...3 5.5 6 11.5 12.5 5.5 4.5
## Percentage...4 4.5 7 11.5 12.5 4.5 3.5
## Percentage...5 5.5 6 11.5 12.5 5.5 4.5
## Percentage...6 5.5 6 11.5 12.5 5.5 4.5
## maxlad_Hdptf maxlad_dptf maxlad_effdist last_Hcbh last_dptf
## Percentage...1 14.5 8 5 6.5 8
## Percentage...2 13.5 8 4 5.5 8
## Percentage...3 11.5 6 4 5.5 6
## Percentage...4 11.5 7 3 4.5 7
## Percentage...5 11.5 6 4 5.5 6
## Percentage...6 11.5 6 4 5.5 6
## last_Hdptf
## Percentage...1 14.5
## Percentage...2 13.5
## Percentage...3 11.5
## Percentage...4 11.5
## Percentage...5 11.5
## Percentage...6 11.5#18. PLOT FUEL LAYERS WITH LAD > 25 % AND CBH BASED ON MAXIMUM LAD PERCENTAGE
library(ggplot2)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)
# Tree metrics derived from get_layers_lad() function
LAD_gt25p_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_gt25perc.txt")
fuels_LAD_metrics <- read.table(LAD_gt25p_path, sep = "\t", header = TRUE)
fuels_LAD_metrics$treeID <- factor(fuels_LAD_metrics$treeID)
trees_name1 <- as.character(fuels_LAD_metrics$treeID)
trees_name2 <- factor(unique(trees_name1))
# Generate plots for fuels LAD metrics
plots_trees_LAD <- get_plots_cbh_LAD(LAD_profiles, fuels_LAD_metrics)
# Save plots for each tree
for (name in names(plots_trees_LAD)) {
plots <- plots_trees_LAD[[name]]
if (!is.null(plots)) {
#print(paste("Saving plot for tree:", name))
ggsave(file.path(system.file("extdata", package = "LadderFuelsR"), paste0( name, "_LAD_25perc", ".tiff")), plot = plots,
width = 6, height = 6, units = "in")
}}
par(mfrow = c(3, 1))
# Plot the first panel
plot(plots_trees_LAD[[1]],width = 6, height = 6)
# Plot the second panel
plot(plots_trees_LAD[[2]],width = 6, height = 6)
# Plot the third panel
plot(plots_trees_LAD[[3]],width = 6, height = 6)
#19. CBH BASED ON THE BREAKING POINT METHOD AND LAD PERCENTAGE
library(dplyr)
library(magrittr)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)
# Tree metrics derived from get_effective_gap() function
distcor_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "6_tree_distances_metrics_corr.txt")
distances_metrics_corr <- read.table(distcor_path, sep = "\t", header = TRUE)
distances_metrics_corr$treeID <- factor(distances_metrics_corr$treeID)
trees_name1 <- as.character(distances_metrics_corr$treeID)
trees_name2 <- factor(unique(trees_name1))
cum_LAD_metrics_list <- list()
for (i in levels(trees_name2)) {
# Filter data for each tree
tree1 <- LAD_profiles |> dplyr::filter(treeID == i)
tree2 <- distances_metrics_corr |> dplyr::filter(treeID == i)
# Get cumulative LAD metrics for each tree
cum_LAD_metrics <- get_cum_break(tree1, tree2)
cum_LAD_metrics_list[[i]] <- cum_LAD_metrics
}## breakpoint estimate(s): 3.118038# Combine the individual data frames
cum_LAD_metrics_all <- dplyr::bind_rows(cum_LAD_metrics_list)
# =======================================================================#
# REORDER COLUMNS
# =======================================================================#
# Get original column names
original_column_names <- colnames(cum_LAD_metrics_all)
# Specify prefixes (adjust accordingly)
prefixes <- c("treeID", "Hcbh", "below", "above", "depth", "Hdepth", "dptf", "Hdptf", "Hdist", "effdist", "max", "last", "cumlad")
# Initialize vector to store new order
new_order <- c()
# Loop over prefixes
for (prefix in prefixes) {
# Find column names matching the current prefix
matching_columns <- grep(paste0("^", prefix), original_column_names, value = TRUE)
# Extract numeric suffixes and order the columns based on these suffixes
numeric_suffixes <- as.numeric(gsub(paste0("^", prefix), "", matching_columns))
matching_columns <- matching_columns[order(numeric_suffixes)]
# Append to new order
new_order <- c(new_order, matching_columns)
}
# Reorder columns
cum_LAD_metrics_order <- cum_LAD_metrics_all[, new_order]
cum_LAD_metrics_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "8_cbh_breaking_point_lad.txt")
write.table(cum_LAD_metrics_order, file = cum_LAD_metrics_path, sep = "\t", row.names = FALSE)
head(cum_LAD_metrics_order)## treeID1 treeID Hcbh1 Hcbh_bp Hcbh1_Hdptf1 below_hcbhbp
## 1 1 1_Eglin_zone1_CROWN 8.487859 8.487859 89.45101 33.925914
## 2 10 10_Eglin_zone1_CROWN 6.809724 6.809724 80.49151 26.699603
## 3 100 100_Eglin_zone1_CROWN 6.843101 6.843101 96.20392 9.172215
## 4 101 101_Eglin_zone1_CROWN NA 7.497583 NA 58.726259
## 5 102 102_Eglin_zone1_CROWN 6.367170 6.367170 81.03784 25.827027
## 6 103 103_Eglin_zone1_CROWN 6.168868 6.168868 92.15937 13.009641
## above_hcbhbp dptf1 Hdptf1 Hdist1 effdist1 maxlad_Hcbh maxlad_Hdptf
## 1 89.45101 7 15.5 7.5 7 8.5 15.5
## 2 80.49151 8 14.5 5.8 6 6.8 14.5
## 3 96.20392 6 12.5 5.8 6 6.8 12.5
## 4 65.31229 NA NA NA NA NA NA
## 5 81.03784 6 12.5 5.4 5 6.4 12.5
## 6 92.15937 6 12.5 5.2 5 6.2 12.5
## maxlad_Hdist maxlad_effdist max_dptf max_Hcbh max_Hdptf max_Hdist max_height
## 1 7.5 7 7 8.5 15.5 7.5 15.5
## 2 5.8 6 8 6.8 14.5 5.8 14.5
## 3 5.8 6 6 6.8 12.5 5.8 12.5
## 4 NA NA NA NA NA NA 12.5
## 5 5.4 5 6 6.4 12.5 5.4 12.5
## 6 5.2 5 6 6.2 12.5 5.2 12.5
## last_effdist last_Hcbh last_Hdptf last_Hdist cumlad
## 1 7 8.5 15.5 7.5 0.4171
## 2 6 6.8 14.5 5.8 0.4013
## 3 6 6.8 12.5 5.8 0.0923
## 4 NA NA NA NA 0.4404
## 5 5 6.4 12.5 5.4 0.2389
## 6 5 6.2 12.5 5.2 0.1120#20. PLOT CBH BASED ON THE BREAKING POINT METHOD AND LAD PERCENTAGE
library(ggplot2)
# LAD profiles derived from normalized ALS data after applying [lad.profile()] function
data_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "LAD_profiles.txt")
LAD_profiles <- read.table(data_path, sep = "\t", header = TRUE)
LAD_profiles$treeID <- factor(LAD_profiles$treeID)
# Tree metrics derived from get_cum_break() function
cum_LAD_metrics_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "8_cbh_breaking_point_lad.txt")
cum_LAD_metrics <- read.table(cum_LAD_metrics_path, sep = "\t", header = TRUE)
cum_LAD_metrics$treeID <- factor(cum_LAD_metrics$treeID)
# Generate cumulative LAD plots
plots_trees_cumlad <- get_plots_cumm(LAD_profiles, cum_LAD_metrics)
# Save plots for each tree
for (name in names(plots_trees_cumlad)) {
plots <- plots_trees_cumlad[[name]]
if (!is.null(plots)) {
#print(paste("Saving plot for tree:", name))
ggsave(file.path(system.file("extdata", package = "LadderFuelsR"), paste0( name, "_cumm_Hcbh", ".tiff")), plot = plots)
}}
par(mfrow = c(3, 1))
# Plot the first panel
plot(plots_trees_cumlad[[1]],width = 6, height = 6)
# Plot the second panel
plot(plots_trees_cumlad[[2]],width = 6, height = 6)
# Plot the third panel
plot(plots_trees_cumlad[[3]],width = 6, height = 6)
#21. JOINING FUEL LADDER PROPERITES WITH CROWN POLYGONS
# Tree metrics derived from get_layers_lad() function
LAD_gt25p_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "7_fuels_lad_gt25perc.txt")
fuels_LAD_metrics <- read.table(LAD_gt25p_path, sep = "\t", header = TRUE)
fuels_LAD_metrics$treeID <- factor(fuels_LAD_metrics$treeID)
# No overlapping crown polygins (output from step 6)
trees_path <- file.path(system.file("extdata", package = "LadderFuelsR"), "crown2m_silva_convex_no_overlap.shp")
tree_file <- st_read(trees_path)## Reading layer `crown2m_silva_convex_no_overlap' from data source
## `D:\OLGA\R_library\LadderFuelsR\extdata\crown2m_silva_convex_no_overlap.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 877 features and 16 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 525912.5 ymin: 3378523 xmax: 526028 ymax: 3378614
## Projected CRS: NAD83 / UTM zone 16Ntree_file$treeID1 <- factor(tree_file$treeID)
crowns_properties<-merge (tree_file,fuels_LAD_metrics, by="treeID1")
# Define a red-to-green color palette
red_to_green_palette <- colorRampPalette(c("red", "green"))
# Plotting with the red-to-green color palette
ggplot() +
geom_sf(data = crowns_properties, aes(fill = Hcbh1)) +
scale_fill_gradientn(colors = red_to_green_palette(100), na.value = "grey50") +
theme_minimal() +
labs(title = "Tree Crowns", fill = "Hcbh1")
crowns_properties$maxlad_Hcbh_factor <- cut(crowns_properties$maxlad_Hcbh, breaks = 5)
# Plotting with a discrete legend
ggplot() +
geom_sf(data = crowns_properties, aes(fill = maxlad_Hcbh_factor)) +
scale_fill_manual(values = red_to_green_palette(5)) +
theme_minimal() +
labs(title = "Tree Crowns", fill = "maxlad_Hcbh")